Cases and social distancing correlation analysis, SCC and SMC

library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(leaflet)
library(lehdr)
library(usmap)
library(lmtest)
library(pracma)
library(lmtest)
library(forecast)
library(vars)
library(rvest)
library(RSelenium)
library(seleniumPipes)
library(dLagM)

options(
  tigris_class = "sf",
  tigris_use_cache = TRUE
)

Get the case data, first for SCC.

# remDr <- remoteDriver(
#   remoteServerAddr = "192.168.86.25",
#   port = 4445L
# )
# remDr$open()
# 
# remDr$navigate("https://app.powerbigov.us/view?r=eyJrIjoiZTg2MTlhMWQtZWE5OC00ZDI3LWE4NjAtMTU3YWYwZDRlOTNmIiwidCI6IjBhYzMyMDJmLWMzZTktNGY1Ni04MzBkLTAxN2QwOWQxNmIzZiJ9")
# 
# webElem <- remDr$findElements(using = "class", value = "column")
# 
# cases <-
#   1:length(webElem) %>%
#   map(function(x){
#     webElem[[x]]$getElementAttribute("aria-label") %>% as.character()
#   }) %>%
#   unlist() %>%
#   as.data.frame()
# 
# scc_cumulative_cases <-
#   cases %>%
#   rename(text = ".") %>%
#   filter(grepl("Total_cases",text)) %>%
#   separate(text, c("date","cases"), sep = "\\.") %>%
#   mutate(
#     date =
#       substr(date,6,nchar(.)) %>%
#       as.Date("%A, %B %d, %Y"),
#     cases =
#       substr(cases,13,nchar(.)) %>%
#       as.numeric()
#   )
# 
# saveRDS(scc_cumulative_cases, file = "/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/scc_cumulative_cases.rds")

scc_cumulative_cases <- readRDS("/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/scc_cumulative_cases.rds")

Also for SMC.

# remDr$navigate("https://app.powerbigov.us/view?r=eyJrIjoiMWI5NmE5M2ItOTUwMC00NGNmLWEzY2UtOTQyODA1YjQ1NWNlIiwidCI6IjBkZmFmNjM1LWEwNGQtNDhjYy1hN2UzLTZkYTFhZjA4ODNmOSJ9")
# 
# webElem <- remDr$findElements(using = "class", value = "column")
# 
# tests <-
#   1:length(webElem) %>%
#   map(function(x){
#     webElem[[x]]$getElementAttribute("aria-label") %>% as.character()
#   }) %>%
#   unlist() %>%
#   as.data.frame()
# 
# tests_clean <-
#   tests %>%
#   rename(text = ".") %>%
#   separate(text, c("date","test_text"), sep = "\\.") %>%
#   separate(test_text, c(NA,"type",NA,"tests")) %>%
#   mutate(
#     date =
#       substr(date,23,nchar(.)) %>%
#       as.Date("%A, %B %d, %Y"),
#     tests =
#       tests %>%
#       as.numeric()
#   ) %>%
#   spread(
#     key = type,
#     value = tests
#   ) %>%
#   mutate(
#     total = Negative + Pending + Positive,
#     perc_positive = Positive/total,
#     perc_positive_movavg = movavg(perc_positive, 7, type = "s")
#   )
# 
# smc_cumulative_cases <- tests_clean %>%
#   mutate(cumulative_cases = cumsum(Positive), cumulative_negative = cumsum(Negative), cumulative_total = cumulative_cases+cumulative_negative, perc_positive_cumulative = cumulative_cases*100 / cumulative_total, perc_positive_cumulative_mov = movavg(perc_positive_cumulative, 7, type = "s"))
# 
# saveRDS(smc_cumulative_cases, file = "/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/smc_cumulative_cases.rds")

smc_cumulative_cases <- readRDS("/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/smc_cumulative_cases.rds")

Get social distancing data.

scc_blockgroups <-
  block_groups("CA","Santa Clara", cb=F, progress_bar=F) %>% 
  st_transform('+proj=longlat +datum=WGS84')

smc_blockgroups <-
  block_groups("CA","San Mateo", cb=F, progress_bar=F) %>% 
  st_transform('+proj=longlat +datum=WGS84')

bay_sd <- readRDS("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Restricted Data Library/Safegraph/covid19analysis/bay_socialdistancing_v2.rds") %>% 
  mutate(date = date_range_start %>%  substr(1,10) %>% as.Date())

# obtaining weekends
weekends <- bay_sd %>% 
  filter(!duplicated(date)) %>% 
  arrange(date) %>% 
  mutate(weekend = ifelse((date %>% as.numeric()) %% 7 %in% c(2,3), T, F)) %>% 
  dplyr::select(date,weekend)

bay_sd <- bay_sd %>% left_join(weekends)

Santa Clara County

scc_cases_sd_daily <- bay_sd %>% 
  filter(origin_census_block_group %in% scc_blockgroups$GEOID) %>%
  group_by(date) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(
    percent_at_home = total_at_home*100/total_devices, 
    percent_leaving_home = (100 - percent_at_home),
  ) %>% 
  left_join(
    scc_cumulative_cases
  ) %>% 
  filter(date >= min(scc_cumulative_cases$date))

# get the baseline percent of people at home
pre_case_growth_at_home_scc <- bay_sd %>%
  filter(date < min(scc_cumulative_cases$date)) %>%
  filter(origin_census_block_group %in% scc_blockgroups$GEOID) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(percent_at_home = total_at_home*100/total_devices, percent_leaving_home = (100 - percent_at_home))

# include change in percent leaving home
scc_cases_sd_daily <- scc_cases_sd_daily %>%
  mutate(leaving_home_dif = percent_leaving_home - pre_case_growth_at_home_scc$percent_leaving_home[1],
         log_cases = log(cases))

# compute number of differences for stationarity
ndiffs(scc_cases_sd_daily$cases)

## [1] 2

ndiffs(scc_cases_sd_daily$log_cases[-1])

## [1] 2

ndiffs(scc_cases_sd_daily$leaving_home_dif)

## [1] 1

scc_test_big <- scc_cases_sd_daily %>%
  mutate(
    dlog_cases = c(NA,diff(log_cases)),
    d2log_cases = c(NA,NA,diff(log_cases, differences = 2)),
    dcases = c(NA,diff(cases)),
    d2cases = c(NA,NA,diff(dcases, differences = 2)),
    dleaving = c(NA,diff(leaving_home_dif)),
    d2leaving = c(NA,NA,diff(leaving_home_dif, differences = 2)),
    cases_mov = movavg(cases, 7, type = "s"),
    log_cases_mov = movavg(log_cases, 7, type = "s"),
    dlog_cases_mov = c(NA,diff(log_cases_mov)),
    d2log_cases_mov = c(NA,NA,diff(log_cases_mov, differences = 2)),
    dcases_mov = c(NA,diff(cases_mov)),
    d2cases_mov = c(NA,diff(dcases_mov)),
    leaving_mov = movavg(leaving_home_dif, 7, type = "s"),
    dleaving_mov = c(NA,diff(leaving_mov)),
    d2leaving_mov = c(NA,diff(dleaving_mov)),
    leaving4 = c(rep(NA,4), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-4)]),
    dleaving4 = c(NA,diff(leaving4)),
    d2leaving4 = c(NA,NA,diff(leaving4, differences = 2)),
    leaving3 = c(rep(NA,3), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-3)]),
    leaving3_mov = movavg(leaving3, 7, type = "s"),
    dleaving3_mov = c(NA,diff(leaving3_mov)),
    d2leaving3_mov = c(NA,NA,diff(leaving3_mov, differences = 2)),
    leaving4_mov = movavg(leaving4, 7, type = "s"),
    dleaving4_mov = c(NA,diff(leaving4_mov)),
    d2leaving4_mov = c(NA,NA,diff(leaving4_mov, differences = 2)),
    leaving5 = c(rep(NA,5), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-5)]),
    leaving5_mov = movavg(leaving5, 7, type = "s"),
    dleaving5_mov = c(NA,diff(leaving5_mov)),
    d2leaving5_mov = c(NA,NA,diff(leaving5_mov, differences = 2)),
    leaving6 = c(rep(NA,6), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-6)]),
    leaving6_mov = movavg(leaving6, 7, type = "s"),
    dleaving6_mov = c(NA,diff(leaving6_mov)),
    d2leaving6_mov = c(NA,NA,diff(leaving6_mov, differences = 2)),
    leaving7 = c(rep(NA,7), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-7)]),
    leaving7_mov = movavg(leaving7, 7, type = "s"),
    dleaving7_mov = c(NA,diff(leaving7_mov)),
    d2leaving7_mov = c(NA,NA,diff(leaving7_mov, differences = 2)),
    leaving8 = c(rep(NA,8), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-8)]),
    leaving8_mov = movavg(leaving8, 7, type = "s"),
    dleaving8_mov = c(NA,diff(leaving8_mov)),
    d2leaving8_mov = c(NA,NA,diff(leaving8_mov, differences = 2)),
    leaving9 = c(rep(NA,9), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-9)]),
    leaving9_mov = movavg(leaving9, 7, type = "s"),
    dleaving9_mov = c(NA,diff(leaving9_mov)),
    d2leaving9_mov = c(NA,NA,diff(leaving9_mov, differences = 2)),
    leaving10 = c(rep(NA,10), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-10)]),
    leaving10_mov = movavg(leaving10, 7, type = "s"),
    dleaving10_mov = c(NA,diff(leaving10_mov)),
    d2leaving10_mov = c(NA,NA,diff(leaving10_mov, differences = 2)),
    past_cases = c(NA, scc_cases_sd_daily$cases[1:(nrow(scc_cases_sd_daily)-1)]), 
    cases_growth_daily = (dcases / past_cases)*100,
    cases_growth_daily_mov = movavg(cases_growth_daily, 7, type = "s")
  ) %>% 
  filter(date >= "2020-03-01")

ndiffs(scc_cases_sd_daily$cases_growth_daily)

## [1] 0

Tests

Granger

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:3) + Lags(dleaving_mov, 1:3)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:3)
##   Res.Df Df      F    Pr(>F)    
## 1     73                        
## 2     76 -3 9.0988 3.424e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:3) + Lags(d2log_cases_mov, 1:3)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     73                 
## 2     76 -3 1.2405 0.3013

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:4) + Lags(dleaving_mov, 1:4)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:4)
##   Res.Df Df      F    Pr(>F)    
## 1     70                        
## 2     74 -4 8.1953 1.764e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:4) + Lags(d2log_cases_mov, 1:4)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     70                 
## 2     74 -4 0.8554 0.4951

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:5) + Lags(dleaving_mov, 1:5)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:5)
##   Res.Df Df      F    Pr(>F)    
## 1     67                        
## 2     72 -5 6.7985 3.506e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:5) + Lags(d2log_cases_mov, 1:5)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:5)
##   Res.Df Df     F Pr(>F)
## 1     67                
## 2     72 -5 0.694 0.6298

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:6) + Lags(dleaving_mov, 1:6)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:6)
##   Res.Df Df      F    Pr(>F)    
## 1     64                        
## 2     70 -6 5.3088 0.0001703 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:6) + Lags(d2log_cases_mov, 1:6)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:6)
##   Res.Df Df     F Pr(>F)
## 1     64                
## 2     70 -6 1.041 0.4074

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:7) + Lags(dleaving_mov, 1:7)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:7)
##   Res.Df Df      F   Pr(>F)   
## 1     61                      
## 2     68 -7 2.9586 0.009776 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:7) + Lags(d2log_cases_mov, 1:7)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:7)
##   Res.Df Df      F  Pr(>F)  
## 1     61                    
## 2     68 -7 1.8237 0.09882 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:8) + Lags(dleaving_mov, 1:8)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:8)
##   Res.Df Df      F    Pr(>F)    
## 1     58                        
## 2     66 -8 4.8424 0.0001338 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:8) + Lags(d2log_cases_mov, 1:8)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:8)
##   Res.Df Df      F  Pr(>F)  
## 1     58                    
## 2     66 -8 2.1953 0.04087 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:9) + Lags(dleaving_mov, 1:9)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:9)
##   Res.Df Df     F    Pr(>F)    
## 1     55                       
## 2     64 -9 3.851 0.0007892 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:9) + Lags(d2log_cases_mov, 1:9)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:9)
##   Res.Df Df      F  Pr(>F)  
## 1     55                    
## 2     64 -9 1.7838 0.09234 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  d2log_cases_mov ~ dleaving_mov,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:10) + Lags(dleaving_mov, 1:10)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:10)
##   Res.Df  Df      F    Pr(>F)    
## 1     52                         
## 2     62 -10 4.0521 0.0003831 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving_mov ~ d2log_cases_mov,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:10) + Lags(d2log_cases_mov, 1:10)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:10)
##   Res.Df  Df      F Pr(>F)
## 1     52                  
## 2     62 -10 1.5036 0.1648

Granger for not moving average

grangertest(
  d2log_cases ~ dleaving,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:3) + Lags(dleaving, 1:3)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:3)
##   Res.Df Df     F   Pr(>F)   
## 1     73                     
## 2     76 -3 4.467 0.006171 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:3) + Lags(d2log_cases, 1:3)
## Model 2: dleaving ~ Lags(dleaving, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     73                 
## 2     76 -3 0.9464 0.4228

grangertest(
  d2log_cases ~ dleaving,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:4) + Lags(dleaving, 1:4)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:4)
##   Res.Df Df      F    Pr(>F)    
## 1     70                        
## 2     74 -4 6.5116 0.0001631 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:4) + Lags(d2log_cases, 1:4)
## Model 2: dleaving ~ Lags(dleaving, 1:4)
##   Res.Df Df    F Pr(>F)
## 1     70               
## 2     74 -4 1.24  0.302

grangertest(
  d2log_cases ~ dleaving,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:5) + Lags(dleaving, 1:5)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:5)
##   Res.Df Df      F   Pr(>F)   
## 1     67                      
## 2     72 -5 4.1674 0.002344 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:5) + Lags(d2log_cases, 1:5)
## Model 2: dleaving ~ Lags(dleaving, 1:5)
##   Res.Df Df      F Pr(>F)
## 1     67                 
## 2     72 -5 0.8606 0.5123

grangertest(
  d2log_cases ~ dleaving,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:6) + Lags(dleaving, 1:6)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:6)
##   Res.Df Df      F  Pr(>F)  
## 1     64                    
## 2     70 -6 2.9787 0.01253 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:6) + Lags(d2log_cases, 1:6)
## Model 2: dleaving ~ Lags(dleaving, 1:6)
##   Res.Df Df     F Pr(>F)
## 1     64                
## 2     70 -6 0.844 0.5408

grangertest(
  d2log_cases ~ dleaving,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:7) + Lags(dleaving, 1:7)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:7)
##   Res.Df Df      F  Pr(>F)  
## 1     61                    
## 2     68 -7 2.7765 0.01423 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:7) + Lags(d2log_cases, 1:7)
## Model 2: dleaving ~ Lags(dleaving, 1:7)
##   Res.Df Df      F Pr(>F)
## 1     61                 
## 2     68 -7 0.9196 0.4978

grangertest(
  d2log_cases ~ dleaving,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:8) + Lags(dleaving, 1:8)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:8)
##   Res.Df Df     F  Pr(>F)  
## 1     58                   
## 2     66 -8 2.765 0.01167 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:8) + Lags(d2log_cases, 1:8)
## Model 2: dleaving ~ Lags(dleaving, 1:8)
##   Res.Df Df      F Pr(>F)
## 1     58                 
## 2     66 -8 1.2574 0.2835

grangertest(
  d2log_cases ~ dleaving,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:9) + Lags(dleaving, 1:9)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:9)
##   Res.Df Df     F  Pr(>F)  
## 1     55                   
## 2     64 -9 2.712 0.01089 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:9) + Lags(d2log_cases, 1:9)
## Model 2: dleaving ~ Lags(dleaving, 1:9)
##   Res.Df Df      F Pr(>F)
## 1     55                 
## 2     64 -9 1.4197 0.2025

grangertest(
  d2log_cases ~ dleaving,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:10) + Lags(dleaving, 1:10)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:10)
##   Res.Df  Df      F  Pr(>F)  
## 1     52                     
## 2     62 -10 2.2658 0.02756 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ d2log_cases,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:10) + Lags(d2log_cases, 1:10)
## Model 2: dleaving ~ Lags(dleaving, 1:10)
##   Res.Df  Df      F Pr(>F)
## 1     52                  
## 2     62 -10 0.9806 0.4715

Granger with case growth rate

grangertest(
  cases_growth_daily ~ dleaving,
  order = 2,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:2) + Lags(dleaving, 1:2)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     76                 
## 2     78 -2 0.8033 0.4516

grangertest(
  dleaving ~ cases_growth_daily,
  order = 2,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:2) + Lags(cases_growth_daily, 1:2)
## Model 2: dleaving ~ Lags(dleaving, 1:2)
##   Res.Df Df      F  Pr(>F)  
## 1     76                    
## 2     78 -2 3.6329 0.03113 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:3) + Lags(dleaving, 1:3)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:3)
##   Res.Df Df      F  Pr(>F)  
## 1     73                    
## 2     76 -3 2.2196 0.09306 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 3,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:3) + Lags(cases_growth_daily, 1:3)
## Model 2: dleaving ~ Lags(dleaving, 1:3)
##   Res.Df Df      F  Pr(>F)  
## 1     73                    
## 2     76 -3 3.6384 0.01663 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:4) + Lags(dleaving, 1:4)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:4)
##   Res.Df Df      F   Pr(>F)   
## 1     70                      
## 2     74 -4 3.7127 0.008478 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 4,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:4) + Lags(cases_growth_daily, 1:4)
## Model 2: dleaving ~ Lags(dleaving, 1:4)
##   Res.Df Df      F  Pr(>F)  
## 1     70                    
## 2     74 -4 2.7896 0.03281 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:5) + Lags(dleaving, 1:5)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:5)
##   Res.Df Df      F  Pr(>F)  
## 1     67                    
## 2     72 -5 3.1335 0.01334 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 5,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:5) + Lags(cases_growth_daily, 1:5)
## Model 2: dleaving ~ Lags(dleaving, 1:5)
##   Res.Df Df      F  Pr(>F)  
## 1     67                    
## 2     72 -5 2.8962 0.01993 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:6) + Lags(dleaving, 1:6)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:6)
##   Res.Df Df      F   Pr(>F)   
## 1     64                      
## 2     70 -6 3.8148 0.002596 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 6,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:6) + Lags(cases_growth_daily, 1:6)
## Model 2: dleaving ~ Lags(dleaving, 1:6)
##   Res.Df Df      F  Pr(>F)  
## 1     64                    
## 2     70 -6 2.6223 0.02457 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:7) + Lags(dleaving, 1:7)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:7)
##   Res.Df Df     F    Pr(>F)    
## 1     61                       
## 2     68 -7 4.484 0.0004429 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 7,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:7) + Lags(cases_growth_daily, 1:7)
## Model 2: dleaving ~ Lags(dleaving, 1:7)
##   Res.Df Df      F Pr(>F)
## 1     61                 
## 2     68 -7 1.7744 0.1089

grangertest(
  cases_growth_daily ~ dleaving,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:8) + Lags(dleaving, 1:8)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:8)
##   Res.Df Df      F  Pr(>F)  
## 1     58                    
## 2     66 -8 2.4973 0.02107 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 8,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:8) + Lags(cases_growth_daily, 1:8)
## Model 2: dleaving ~ Lags(dleaving, 1:8)
##   Res.Df Df      F  Pr(>F)  
## 1     58                    
## 2     66 -8 2.5929 0.01706 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:9) + Lags(dleaving, 1:9)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:9)
##   Res.Df Df      F  Pr(>F)  
## 1     55                    
## 2     64 -9 1.8937 0.07217 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 9,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:9) + Lags(cases_growth_daily, 1:9)
## Model 2: dleaving ~ Lags(dleaving, 1:9)
##   Res.Df Df      F Pr(>F)  
## 1     55                   
## 2     64 -9 2.2668 0.0307 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  cases_growth_daily ~ dleaving,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:10) + Lags(dleaving, 1:10)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:10)
##   Res.Df  Df      F  Pr(>F)  
## 1     52                     
## 2     62 -10 1.9322 0.06139 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grangertest(
  dleaving ~ cases_growth_daily,
  order = 10,
  data = scc_test_big
)

## Granger causality test
## 
## Model 1: dleaving ~ Lags(dleaving, 1:10) + Lags(cases_growth_daily, 1:10)
## Model 2: dleaving ~ Lags(dleaving, 1:10)
##   Res.Df  Df      F   Pr(>F)   
## 1     52                       
## 2     62 -10 3.5695 0.001183 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Also try ardlDlm

# 4th order lag on x only
testing_ardl4 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c()))
summary(testing_ardl4)

## 
## Time series regression with "ts" data:
## Start = 5, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0213610 -0.0014294  0.0006449  0.0016122  0.0266611 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0009501  0.0008704  -1.092    0.279    
## X.4         -0.0006286  0.0012117  -0.519    0.605    
## Y.1         -0.0364700  0.1170492  -0.312    0.756    
## Y.2          0.2251433  0.1197267   1.880    0.064 .  
## Y.3          0.4227324  0.1001782   4.220 6.94e-05 ***
## Y.4          0.0537929  0.1099349   0.489    0.626    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007004 on 73 degrees of freedom
## Multiple R-squared:  0.2386, Adjusted R-squared:  0.1865 
## F-statistic: 4.576 on 5 and 73 DF,  p-value: 0.001098

testing_ardl4 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c()))
summary(testing_ardl4)

## 
## Time series regression with "ts" data:
## Start = 5, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.110146 -0.008604  0.001436  0.009713  0.095544 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.005328   0.003753  -1.420  0.15996    
## X.4          0.001879   0.001196   1.572  0.12026    
## Y.1         -0.648722   0.114527  -5.664 2.75e-07 ***
## Y.2         -0.575380   0.123212  -4.670 1.34e-05 ***
## Y.3         -0.419613   0.128582  -3.263  0.00168 ** 
## Y.4         -0.166485   0.093575  -1.779  0.07938 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03261 on 73 degrees of freedom
## Multiple R-squared:  0.3428, Adjusted R-squared:  0.2978 
## F-statistic: 7.616 on 5 and 73 DF,  p-value: 8.376e-06

# 5th order lag on x only
testing_ardl5 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 5, q = 5, remove = list(p = c(0,1,2,3,4), q=c()))
summary(testing_ardl5)

## 
## Time series regression with "ts" data:
## Start = 6, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0222745 -0.0009438  0.0008675  0.0020159  0.0155077 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.0014185  0.0007383  -1.921   0.0587 .
## X.5          0.0011167  0.0010129   1.102   0.2740  
## Y.1          0.1180938  0.0976595   1.209   0.2306  
## Y.2         -0.0262672  0.0977324  -0.269   0.7889  
## Y.3          0.1833054  0.1023435   1.791   0.0775 .
## Y.4         -0.0544539  0.0932886  -0.584   0.5613  
## Y.5          0.1538378  0.0919230   1.674   0.0986 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005844 on 71 degrees of freedom
## Multiple R-squared:  0.2397, Adjusted R-squared:  0.1754 
## F-statistic:  3.73 on 6 and 71 DF,  p-value: 0.002788

testing_ardl5 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 5, q = 5, remove = list(p = c(0,1,2,3,4), q=c()))
summary(testing_ardl5)

## 
## Time series regression with "ts" data:
## Start = 6, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.101496 -0.006203  0.002880  0.011948  0.080923 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.006203   0.003331  -1.862   0.0667 .  
## X.5         -0.002759   0.001062  -2.598   0.0114 *  
## Y.1         -0.582852   0.101810  -5.725  2.3e-07 ***
## Y.2         -0.312538   0.119545  -2.614   0.0109 *  
## Y.3         -0.311636   0.122142  -2.551   0.0129 *  
## Y.4          0.199758   0.119726   1.668   0.0996 .  
## Y.5          0.214992   0.083142   2.586   0.0118 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02836 on 71 degrees of freedom
## Multiple R-squared:  0.4913, Adjusted R-squared:  0.4483 
## F-statistic: 11.43 on 6 and 71 DF,  p-value: 6.638e-09

# 6th order lag on x only
testing_ardl6 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 6, q = 6, remove = list(p = c(0,1,2,3,4,5), q=c()))
summary(testing_ardl6)

## 
## Time series regression with "ts" data:
## Start = 7, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0208017 -0.0008285  0.0008797  0.0017561  0.0147256 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.0014034  0.0007658  -1.833   0.0712 .
## X.6          0.0012260  0.0010361   1.183   0.2407  
## Y.1          0.1518539  0.1189929   1.276   0.2062  
## Y.2          0.0062912  0.0988165   0.064   0.9494  
## Y.3          0.2571029  0.0983026   2.615   0.0109 *
## Y.4         -0.0797924  0.1055060  -0.756   0.4521  
## Y.5          0.1468134  0.0940073   1.562   0.1229  
## Y.6         -0.1412984  0.0939433  -1.504   0.1371  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005852 on 69 degrees of freedom
## Multiple R-squared:  0.2486, Adjusted R-squared:  0.1724 
## F-statistic: 3.261 on 7 and 69 DF,  p-value: 0.004746

testing_ardl6 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 6, q = 6, remove = list(p = c(0,1,2,3,4,5), q=c()))
summary(testing_ardl6)

## 
## Time series regression with "ts" data:
## Start = 7, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.084402 -0.003895  0.005013  0.010127  0.065618 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0091967  0.0029874  -3.078 0.002984 ** 
## X.6          0.0004921  0.0010058   0.489 0.626196    
## Y.1         -0.5913184  0.1037030  -5.702 2.70e-07 ***
## Y.2         -0.4276185  0.1072323  -3.988 0.000164 ***
## Y.3         -0.5602566  0.1091638  -5.132 2.52e-06 ***
## Y.4         -0.1309827  0.1114987  -1.175 0.244135    
## Y.5         -0.3088103  0.1074100  -2.875 0.005367 ** 
## Y.6         -0.3795404  0.0769707  -4.931 5.42e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02471 on 69 degrees of freedom
## Multiple R-squared:  0.6223, Adjusted R-squared:  0.584 
## F-statistic: 16.24 on 7 and 69 DF,  p-value: 1.953e-12

# 7th order lag on x only
testing_ardl7 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 7, q = 7, remove = list(p = c(0,1,2,3,4,5,6), q=c()))
summary(testing_ardl7)

## 
## Time series regression with "ts" data:
## Start = 8, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0229043 -0.0012596  0.0003516  0.0013377  0.0111903 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.0006688  0.0006501  -1.029  0.30727   
## X.7          0.0009612  0.0008592   1.119  0.26725   
## Y.1          0.0922179  0.0988767   0.933  0.35435   
## Y.2          0.3199778  0.0988280   3.238  0.00188 **
## Y.3          0.1504645  0.0811417   1.854  0.06809 . 
## Y.4          0.1045992  0.0846600   1.236  0.22095   
## Y.5          0.1618351  0.0870405   1.859  0.06738 . 
## Y.6         -0.0780031  0.0785256  -0.993  0.32412   
## Y.7         -0.2322415  0.0783790  -2.963  0.00421 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004804 on 67 degrees of freedom
## Multiple R-squared:  0.4223, Adjusted R-squared:  0.3534 
## F-statistic: 6.123 on 8 and 67 DF,  p-value: 6.573e-06

testing_ardl7 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 7, q = 7, remove = list(p = c(0,1,2,3,4,5,6), q=c()))
summary(testing_ardl7)

## 
## Time series regression with "ts" data:
## Start = 8, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.057916 -0.005367  0.001760  0.007000  0.048366 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0031241  0.0024762  -1.262  0.21146    
## X.7          0.0019605  0.0007771   2.523  0.01402 *  
## Y.1         -0.3011860  0.0928512  -3.244  0.00184 ** 
## Y.2         -0.1065801  0.0969607  -1.099  0.27561    
## Y.3         -0.3201834  0.0917558  -3.490  0.00086 ***
## Y.4          0.1257686  0.0991218   1.269  0.20889    
## Y.5         -0.0768311  0.0870402  -0.883  0.38055    
## Y.6         -0.1397452  0.0878037  -1.592  0.11619    
## Y.7          0.2323405  0.0691169   3.362  0.00128 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01905 on 67 degrees of freedom
## Multiple R-squared:  0.7695, Adjusted R-squared:  0.742 
## F-statistic: 27.96 on 8 and 67 DF,  p-value: < 2.2e-16

# all
testing_ardl = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 15, q = 15)
summary(testing_ardl)

## 
## Time series regression with "ts" data:
## Start = 16, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0047279 -0.0010529 -0.0001463  0.0012977  0.0040302 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.0012875  0.0004644  -2.773  0.00875 **
## X.t          0.0017693  0.0008860   1.997  0.05343 . 
## X.1          0.0009649  0.0010808   0.893  0.37791   
## X.2          0.0005149  0.0009821   0.524  0.60328   
## X.3         -0.0005939  0.0010123  -0.587  0.56105   
## X.4          0.0002089  0.0009913   0.211  0.83429   
## X.5         -0.0002293  0.0009836  -0.233  0.81701   
## X.6         -0.0001107  0.0010769  -0.103  0.91872   
## X.7          0.0012352  0.0011448   1.079  0.28777   
## X.8          0.0001928  0.0012221   0.158  0.87555   
## X.9          0.0002823  0.0010341   0.273  0.78640   
## X.10         0.0020642  0.0010017   2.061  0.04662 * 
## X.11        -0.0002756  0.0010028  -0.275  0.78499   
## X.12        -0.0002614  0.0009913  -0.264  0.79356   
## X.13        -0.0021029  0.0009363  -2.246  0.03093 * 
## X.14         0.0028202  0.0009218   3.059  0.00417 **
## X.15         0.0001016  0.0009576   0.106  0.91608   
## Y.1          0.1151045  0.1442918   0.798  0.43026   
## Y.2          0.1030344  0.1355770   0.760  0.45222   
## Y.3          0.0148629  0.1346612   0.110  0.91273   
## Y.4          0.1182629  0.1258654   0.940  0.35369   
## Y.5         -0.1212860  0.1159990  -1.046  0.30273   
## Y.6         -0.1238151  0.1064913  -1.163  0.25261   
## Y.7         -0.2529957  0.1142312  -2.215  0.03319 * 
## Y.8         -0.0526962  0.0976121  -0.540  0.59262   
## Y.9          0.0865780  0.0809040   1.070  0.29168   
## Y.10        -0.0811094  0.0794617  -1.021  0.31419   
## Y.11        -0.0081931  0.0665314  -0.123  0.90268   
## Y.12         0.0450652  0.0639009   0.705  0.48520   
## Y.13         0.0159947  0.0577512   0.277  0.78340   
## Y.14        -0.0224993  0.0547341  -0.411  0.68346   
## Y.15        -0.1010432  0.0603504  -1.674  0.10274   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.00222 on 36 degrees of freedom
## Multiple R-squared:  0.9006, Adjusted R-squared:  0.815 
## F-statistic: 10.52 on 31 and 36 DF,  p-value: 1.511e-10

testing_ardl = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 15, q = 15)
summary(testing_ardl)

## 
## Time series regression with "ts" data:
## Start = 16, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.016155 -0.004704 -0.001154  0.006352  0.015003 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.293e-03  1.941e-03  -1.697  0.09841 .  
## X.t          2.027e-03  6.505e-04   3.116  0.00359 ** 
## X.1          1.500e-03  6.940e-04   2.162  0.03736 *  
## X.2          2.380e-03  6.940e-04   3.430  0.00153 ** 
## X.3          5.943e-04  7.676e-04   0.774  0.44383    
## X.4          1.412e-03  7.068e-04   1.997  0.05338 .  
## X.5          1.357e-03  7.601e-04   1.786  0.08261 .  
## X.6          1.555e-04  6.938e-04   0.224  0.82397    
## X.7          9.477e-04  6.627e-04   1.430  0.16132    
## X.8          7.853e-05  6.578e-04   0.119  0.90563    
## X.9          1.075e-03  6.735e-04   1.596  0.11919    
## X.10         2.143e-03  6.861e-04   3.123  0.00352 ** 
## X.11         1.606e-03  7.094e-04   2.264  0.02972 *  
## X.12         1.580e-03  7.038e-04   2.244  0.03104 *  
## X.13        -2.323e-04  7.488e-04  -0.310  0.75823    
## X.14         2.022e-03  7.047e-04   2.869  0.00684 ** 
## X.15         3.584e-04  7.411e-04   0.484  0.63160    
## Y.1         -6.005e-01  1.615e-01  -3.718  0.00068 ***
## Y.2         -4.422e-01  1.447e-01  -3.056  0.00421 ** 
## Y.3         -3.223e-01  1.518e-01  -2.123  0.04070 *  
## Y.4         -2.256e-01  1.451e-01  -1.554  0.12883    
## Y.5         -1.919e-01  1.332e-01  -1.441  0.15830    
## Y.6         -6.598e-02  1.229e-01  -0.537  0.59463    
## Y.7         -7.444e-02  1.116e-01  -0.667  0.50889    
## Y.8         -1.354e-01  1.057e-01  -1.281  0.20825    
## Y.9         -1.528e-01  9.428e-02  -1.621  0.11380    
## Y.10        -2.657e-01  9.598e-02  -2.769  0.00884 ** 
## Y.11        -2.394e-01  1.028e-01  -2.330  0.02554 *  
## Y.12        -9.125e-02  1.042e-01  -0.876  0.38695    
## Y.13         5.172e-02  9.874e-02   0.524  0.60359    
## Y.14         1.025e-01  8.732e-02   1.174  0.24804    
## Y.15         1.121e-01  6.369e-02   1.760  0.08683 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0101 on 36 degrees of freedom
## Multiple R-squared:  0.7763, Adjusted R-squared:  0.5837 
## F-statistic:  4.03 on 31 and 36 DF,  p-value: 4.468e-05

# try with not log of cases
testing_ardl_not_log = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2cases, p = 15, q = 15)
summary(testing_ardl_not_log)

## 
## Time series regression with "ts" data:
## Start = 16, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.3602  -7.9741  -0.3545   5.9212  17.7407 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.75083    1.63037   0.461 0.647909    
## X.t          0.67083    0.62922   1.066 0.293465    
## X.1          0.34348    0.60532   0.567 0.573945    
## X.2          0.83229    0.55930   1.488 0.145437    
## X.3         -0.42042    0.60086  -0.700 0.488610    
## X.4          0.75119    0.59930   1.253 0.218126    
## X.5         -0.62798    0.62731  -1.001 0.323479    
## X.6         -0.20905    0.60994  -0.343 0.733792    
## X.7         -0.80352    0.59821  -1.343 0.187608    
## X.8         -0.30680    0.60241  -0.509 0.613657    
## X.9          0.07638    0.58368   0.131 0.896615    
## X.10         1.22742    0.57486   2.135 0.039623 *  
## X.11         0.47862    0.59857   0.800 0.429184    
## X.12         0.26747    0.59655   0.448 0.656584    
## X.13        -0.99414    0.58608  -1.696 0.098471 .  
## X.14         3.04977    0.60491   5.042 1.32e-05 ***
## X.15         0.27318    0.74302   0.368 0.715273    
## Y.1         -1.31077    0.16042  -8.171 1.02e-09 ***
## Y.2         -1.47673    0.23695  -6.232 3.40e-07 ***
## Y.3         -1.76798    0.29986  -5.896 9.57e-07 ***
## Y.4         -2.03421    0.36643  -5.551 2.77e-06 ***
## Y.5         -1.86483    0.42127  -4.427 8.53e-05 ***
## Y.6         -1.71959    0.44075  -3.901 0.000402 ***
## Y.7         -1.46637    0.45636  -3.213 0.002767 ** 
## Y.8         -1.29391    0.45227  -2.861 0.006992 ** 
## Y.9         -1.12452    0.43699  -2.573 0.014337 *  
## Y.10        -1.00497    0.40279  -2.495 0.017320 *  
## Y.11        -0.86689    0.36373  -2.383 0.022553 *  
## Y.12        -0.74504    0.32034  -2.326 0.025778 *  
## Y.13        -0.56437    0.27311  -2.066 0.046038 *  
## Y.14        -0.23432    0.20662  -1.134 0.264275    
## Y.15        -0.11322    0.12419  -0.912 0.368015    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.48 on 36 degrees of freedom
## Multiple R-squared:  0.9013, Adjusted R-squared:  0.8162 
## F-statistic:  10.6 on 31 and 36 DF,  p-value: 1.35e-10

# pre april 15th
scc_test_pre_april15 <- scc_test_big %>%
  filter(date <=  "2020-04-15")

# all with this dataset
testing_ardl_pre = ardlDlm(x = scc_test_pre_april15$dleaving_mov, y = scc_test_pre_april15$d2log_cases_mov, p = 4, q = 4, remove = list(p = c(0,1,2,3)))
summary(testing_ardl_pre)

## 
## Time series regression with "ts" data:
## Start = 5, End = 46
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0193197 -0.0028480 -0.0002416  0.0022492  0.0268557 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.002229   0.001805  -1.235  0.22479   
## X.4         -0.001360   0.001975  -0.689  0.49539   
## Y.1         -0.044068   0.168114  -0.262  0.79471   
## Y.2          0.241224   0.176943   1.363  0.18126   
## Y.3          0.430497   0.143925   2.991  0.00499 **
## Y.4          0.074366   0.157262   0.473  0.63915   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009797 on 36 degrees of freedom
## Multiple R-squared:  0.2273, Adjusted R-squared:   0.12 
## F-statistic: 2.118 on 5 and 36 DF,  p-value: 0.08554

testing_ardl_pre = ardlDlm(x = scc_test_pre_april15$dleaving, y = scc_test_pre_april15$d2log_cases, p = 4, q = 4, remove = list(p = c(0,1,2,3)))
summary(testing_ardl_pre)

## 
## Time series regression with "ts" data:
## Start = 5, End = 46
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.101052 -0.023588 -0.003848  0.018214  0.096246 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.009740   0.007306  -1.333 0.190892    
## X.4          0.002338   0.001942   1.204 0.236427    
## Y.1         -0.673451   0.162609  -4.142 0.000199 ***
## Y.2         -0.614642   0.176500  -3.482 0.001322 ** 
## Y.3         -0.461037   0.185506  -2.485 0.017726 *  
## Y.4         -0.197461   0.137284  -1.438 0.158975    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04549 on 36 degrees of freedom
## Multiple R-squared:  0.3602, Adjusted R-squared:  0.2713 
## F-statistic: 4.053 on 5 and 36 DF,  p-value: 0.005082

# ardlBoundOrders
ardlBoundOrders(data = scc_test_big %>% dplyr::select(dleaving, d2log_cases), d2log_cases~dleaving)

## $p
##   dleaving
## 1       13
## 
## $q
## [1] 15
## 
## $Stat.table
##            q = 1     q = 2     q = 3     q = 4     q = 5     q = 6     q = 7
## p = 1  -311.9603 -335.3210 -330.1968 -345.8278 -366.9638 -385.0672 -376.6914
## p = 2  -320.4750 -335.3670 -332.4970 -345.0491 -367.8052 -386.9875 -378.9446
## p = 3  -330.9842 -330.9842 -337.2451 -345.8034 -367.3441 -385.7307 -377.8713
## p = 4  -322.4556 -343.6761 -343.6761 -345.1024 -365.3940 -383.7681 -375.8753
## p = 5  -321.1203 -348.4434 -346.5467 -346.5467 -363.7819 -381.8387 -374.0109
## p = 6  -347.5903 -371.9278 -369.9330 -374.0278 -374.0278 -382.2868 -374.4121
## p = 7  -340.6690 -377.1889 -376.5889 -374.8685 -375.0555 -375.0555 -373.6369
## p = 8  -354.3639 -366.6705 -364.7182 -363.3313 -374.9599 -375.8398 -375.8398
## p = 9  -350.2095 -364.4795 -362.8167 -364.8244 -372.7247 -374.0654 -379.9239
## p = 10 -369.6317 -374.5996 -372.9301 -378.0771 -379.3319 -377.5179 -377.8421
## p = 11 -345.7944 -357.9215 -356.0887 -355.1563 -364.9176 -363.8470 -365.3090
## p = 12 -358.2348 -365.2979 -365.4377 -373.5496 -373.7836 -372.2457 -371.6928
## p = 13 -392.7529 -396.8967 -402.1087 -404.1200 -404.5720 -406.0057 -405.0748
## p = 14 -388.2767 -395.8255 -398.5145 -406.6047 -404.6167 -402.7657 -400.9487
## p = 15 -401.6321 -402.1071 -410.1066 -412.1899 -411.2751 -411.8378 -411.9559
##            q = 8     q = 9    q = 10    q = 11    q = 12    q = 13    q = 14
## p = 1  -384.8290 -378.0385 -373.1607 -373.3674 -376.3853 -374.1426 -382.4685
## p = 2  -386.4006 -379.9632 -373.8916 -372.9623 -375.3222 -372.7886 -380.6993
## p = 3  -387.5460 -380.4194 -375.1833 -373.8577 -378.5349 -376.0810 -381.5210
## p = 4  -387.6737 -381.6816 -376.2598 -373.6625 -377.9804 -375.6735 -383.3153
## p = 5  -387.7739 -381.9738 -375.2375 -373.2163 -376.4734 -373.7979 -383.0741
## p = 6  -386.4163 -380.7912 -374.2681 -371.6942 -375.4924 -373.1138 -384.3340
## p = 7  -384.4275 -378.8401 -372.2728 -370.3987 -373.5037 -371.2976 -382.3836
## p = 8  -385.9502 -379.2483 -373.1510 -369.4446 -371.6045 -369.5184 -380.5363
## p = 9  -379.9239 -381.0504 -375.6413 -373.6736 -373.2675 -374.9434 -383.6552
## p = 10 -377.6418 -377.6418 -377.6468 -379.1690 -376.8519 -377.2439 -389.2916
## p = 11 -369.6717 -371.1493 -371.1493 -382.4594 -381.5857 -381.1412 -397.1622
## p = 12 -371.6109 -370.0949 -381.8265 -381.8265 -379.8830 -383.2483 -399.0410
## p = 13 -403.1759 -404.1079 -406.4798 -410.7621 -410.7621 -409.4577 -410.7312
## p = 14 -399.1349 -399.6236 -401.1789 -407.4241 -405.5559 -405.5559 -409.1715
## p = 15 -409.9693 -413.1741 -414.0036 -417.9835 -418.3189 -420.3182 -420.3182
##           q = 15
## p = 1  -407.8020
## p = 2  -405.8318
## p = 3  -404.0462
## p = 4  -402.3661
## p = 5  -402.8153
## p = 6  -401.2525
## p = 7  -399.9218
## p = 8  -397.9698
## p = 9  -398.4864
## p = 10 -403.0634
## p = 11 -406.7364
## p = 12 -407.4048
## p = 13 -424.1335
## p = 14 -423.2693
## p = 15 -421.5154
## 
## $min.Stat
## [1] -424.1335

ARDL test might not require stationarity?

# 4th order lag on x only, original variables
testing_ardl4_2 = ardlDlm(x = scc_test_big$leaving_home_dif, y = scc_test_big$cases, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c(1,2,3,4)))
summary(testing_ardl4_2)

## 
## Time series regression with "ts" data:
## Start = 5, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -7.302e-14 -6.010e-16  2.620e-16  1.336e-15  2.169e-14 
## 
## Coefficients:
##               Estimate Std. Error    t value Pr(>|t|)    
## (Intercept) -1.117e-13  3.193e-15 -3.497e+01  < 2e-16 ***
## X.4          3.170e-15  1.988e-16  1.595e+01  < 2e-16 ***
## Y.1         -5.704e-16  6.813e-17 -8.373e+00 2.33e-12 ***
## Y.0          1.000e+00  6.988e-17  1.431e+16  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.146e-15 on 75 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.242e+35 on 3 and 75 DF,  p-value: < 2.2e-16

# ardlBoundOrders
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, cases), cases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               15
## 
## $q
## [1] 5
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  650.2973 641.6270 635.7257 630.1846 625.1821 619.8513 607.7379 595.1162
## p = 2  641.7394 643.1781 636.6205 630.6678 625.6807 620.4806 609.3495 596.4008
## p = 3  636.3960 636.3960 638.3134 632.5898 627.6447 622.4095 611.1910 598.4001
## p = 4  629.9057 631.4029 631.4029 632.5275 626.8448 621.1701 610.8961 596.5140
## p = 5  623.7879 624.9465 626.8140 626.8140 628.8048 623.1501 612.7791 597.5042
## p = 6  615.5902 617.0908 618.9111 620.3378 620.3378 621.4988 612.8007 598.1945
## p = 7  610.4319 612.2153 613.6291 614.5274 616.3115 616.3115 614.6655 600.1745
## p = 8  604.3551 606.2545 606.9379 606.6977 607.7826 608.7765 608.7765 601.2606
## p = 9  593.0528 594.7411 595.0907 596.0349 598.0311 599.9950 590.3761 590.3761
## p = 10 587.3863 588.8949 589.3714 589.8791 591.8627 593.6131 585.9210 586.9302
## p = 11 580.3251 581.4935 582.6880 583.9051 585.6084 587.2188 579.6786 580.3983
## p = 12 574.7763 575.4343 576.9362 578.4911 580.0743 581.3619 575.0197 575.4812
## p = 13 560.2912 561.3678 559.4739 557.5347 558.9601 560.6631 560.3960 559.7173
## p = 14 534.3470 535.9078 533.5057 535.0944 529.8075 531.7172 531.1798 532.0997
## p = 15 527.2519 529.1652 525.2829 526.4698 522.1444 523.8479 523.4097 523.3561
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  588.1154 583.1995 576.7651 571.5399 566.1110 559.6037 550.7439
## p = 2  589.1846 583.7996 576.8848 571.8229 566.3383 559.7792 552.2000
## p = 3  591.1773 585.7953 578.7397 573.6795 568.2615 561.7113 554.1987
## p = 4  589.3260 584.0780 577.8451 573.0150 567.6606 560.8485 553.7741
## p = 5  589.4609 583.8559 578.1385 573.5110 568.4760 561.5495 554.1507
## p = 6  591.0250 585.6184 579.8033 575.1108 569.9874 563.0697 555.6178
## p = 7  593.0147 587.5008 581.5781 576.9497 571.8665 565.0507 557.5126
## p = 8  593.1013 587.2463 580.8001 576.2845 570.9630 564.8165 557.0980
## p = 9  592.0671 586.7181 580.8421 576.1842 570.6466 563.5534 555.0586
## p = 10 586.9302 588.3415 582.7671 578.0914 572.3785 564.7322 556.1961
## p = 11 582.3082 582.3082 584.1915 579.5716 574.1144 566.7244 558.0308
## p = 12 577.4601 579.3827 579.3827 581.3374 576.0303 568.7184 559.6443
## p = 13 561.5990 562.8479 564.6884 564.6884 564.2771 552.9216 547.1332
## p = 14 534.0976 535.3977 537.3927 538.9529 538.9529 539.5749 527.2618
## p = 15 525.2618 525.4173 527.2359 527.3845 529.2833 529.2833 529.0907
## 
## $min.Stat
## [1] 522.1444

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, dcases), dcases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               15
## 
## $q
## [1] 6
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  665.0170 658.5016 651.9311 644.0689 637.0321 625.7158 612.0242 606.5196
## p = 2  658.2912 659.8454 653.6515 645.8323 638.7672 626.9381 613.5689 608.0718
## p = 3  650.2567 650.2567 651.8360 644.9846 637.5948 625.4092 613.6611 607.9474
## p = 4  645.1425 646.8280 646.8280 646.9725 639.5079 627.4010 615.5358 609.7576
## p = 5  635.8497 637.7094 639.7091 639.7091 639.9034 628.1471 616.9887 611.4031
## p = 6  621.2309 623.2291 625.1908 626.1594 626.1594 622.8506 611.5754 606.7288
## p = 7  610.4581 612.3972 614.3961 615.4649 616.7337 616.7337 611.6690 606.9245
## p = 8  618.4858 619.9157 620.8035 618.4498 615.7614 606.9768 606.9768 608.8849
## p = 9  605.2032 607.2027 608.0947 605.9398 605.3075 597.4828 595.7400 595.7400
## p = 10 599.3769 601.1624 602.6453 599.8818 598.1908 589.9821 589.3357 591.2508
## p = 11 593.8135 595.8108 597.6681 594.5465 593.2131 585.4305 584.2312 586.2021
## p = 12 583.0040 584.8137 586.3662 586.5914 586.4300 578.9500 577.8757 579.8754
## p = 13 559.2451 558.7821 560.3956 555.6145 557.1194 554.9877 555.6511 557.3983
## p = 14 536.6600 537.0920 538.6950 525.3902 527.0979 527.2312 528.4847 530.4034
## p = 15 537.3860 538.1000 539.0681 524.9590 526.7615 523.3487 524.8220 526.4190
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  600.8366 595.7055 589.1558 584.2013 578.7239 570.7727 562.7727
## p = 2  602.3376 597.3340 590.9455 585.9195 580.3933 572.5229 564.3235
## p = 3  601.8537 596.7337 591.1873 586.1220 580.1873 572.5213 563.7332
## p = 4  603.7968 598.7284 593.1797 588.0903 582.1448 574.5212 565.7221
## p = 5  605.5864 600.4345 594.5913 589.5609 583.3498 575.6893 567.6035
## p = 6  601.5609 596.3216 589.6206 584.3496 578.2406 570.8847 562.5516
## p = 7  601.9416 596.6944 590.0650 584.3657 577.7259 570.8763 562.9464
## p = 8  603.9337 598.6662 592.0493 586.3156 579.5229 572.8586 564.9062
## p = 9  597.7292 592.4713 585.5721 579.7418 571.0319 565.3131 556.3717
## p = 10 591.2508 592.9974 586.6565 580.4958 570.7486 565.6624 556.6824
## p = 11 588.1531 588.1531 588.6562 582.4936 572.5863 567.6060 558.6622
## p = 12 581.7917 583.7872 583.7872 584.2480 574.3161 569.4655 560.6430
## p = 13 559.3781 560.8696 562.2773 562.2773 558.5230 553.2529 547.6419
## p = 14 532.2670 534.0656 536.0643 537.6763 537.6763 537.5414 528.6798
## p = 15 528.3661 529.4351 531.3052 530.2919 528.3901 528.3901 529.7897
## 
## $min.Stat
## [1] 523.3487

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, d2cases), d2cases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               14
## 
## $q
## [1] 15
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  708.0255 695.0573 683.2586 675.0720 638.6222 630.4959 623.7866 618.6388
## p = 2  700.9680 696.4089 684.8028 676.9067 640.6171 632.4896 625.7866 620.6376
## p = 3  687.3660 687.3660 684.3979 676.4178 641.8226 633.8255 627.3596 622.1676
## p = 4  686.2383 678.4874 678.4874 678.0445 643.2502 634.9395 628.3330 623.0245
## p = 5  680.7179 676.9110 672.2588 672.2588 645.1659 636.8969 630.3307 625.0236
## p = 6  669.9032 666.0111 664.7156 636.8153 636.8153 635.7741 629.1239 624.0158
## p = 7  657.4764 654.9395 651.9999 648.5933 630.6631 630.6631 631.1187 626.0006
## p = 8  659.8595 656.9329 653.5996 652.5340 626.0472 625.0618 625.0618 626.8475
## p = 9  652.4397 649.8297 646.0875 645.4320 619.9217 619.8871 620.3688 620.3688
## p = 10 645.4900 643.3666 639.8260 639.3664 613.5580 613.7928 613.8767 615.5973
## p = 11 637.2523 634.8880 631.7602 630.9255 605.3217 604.6739 604.8761 606.5981
## p = 12 629.9253 628.4757 625.3575 625.1048 599.8283 599.1141 598.7457 600.2866
## p = 13 612.6199 612.3519 609.7855 604.1542 587.2176 587.0744 587.4587 589.1213
## p = 14 583.0053 584.3845 568.4426 562.1544 552.3551 549.5834 546.9471 547.9067
## p = 15 574.3285 575.4507 557.0499 558.0778 545.1261 542.9300 541.0397 542.3738
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  613.0900 607.5343 602.2246 589.6840 581.7467 576.0985 569.5043
## p = 2  615.0836 609.5210 604.2234 591.6387 583.7196 578.0765 571.4873
## p = 3  616.4500 610.8603 605.7660 592.8897 584.6454 579.3115 572.5631
## p = 4  617.5113 611.9342 606.9250 594.4230 586.3950 581.1269 574.0788
## p = 5  619.4889 613.9119 608.9097 596.4163 588.3581 583.0644 576.0767
## p = 6  618.7869 613.1851 608.1755 595.9841 587.4106 581.8646 574.6438
## p = 7  620.7803 615.1715 610.1700 597.9841 589.4096 583.8613 576.5998
## p = 8  621.6175 616.0894 611.0284 598.5636 589.6877 583.9511 577.2529
## p = 9  622.1537 616.1944 611.2535 598.6515 589.8647 584.4256 577.3277
## p = 10 615.5973 617.2211 612.0947 599.4522 590.7582 585.4951 578.3009
## p = 11 608.5068 608.5068 610.5032 597.8410 589.1579 583.9348 575.2132
## p = 12 602.2501 604.2497 604.2497 598.2797 589.5414 584.0929 573.7731
## p = 13 591.1199 593.0929 588.6314 588.6314 587.5166 582.2531 573.8226
## p = 14 548.2993 549.8267 550.7611 547.5938 547.5938 547.9128 531.2568
## p = 15 542.3616 544.2568 546.0354 540.2210 539.2820 539.2820 533.1557
## 
## $min.Stat
## [1] 531.2568

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, d2log_cases), d2log_cases~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               14
## 
## $q
## [1] 15
## 
## $Stat.table
##            q = 1     q = 2     q = 3     q = 4     q = 5     q = 6     q = 7
## p = 1  -308.6309 -330.6783 -326.9035 -342.4762 -363.8391 -385.6209 -377.3186
## p = 2  -323.3534 -338.2787 -337.0203 -345.5516 -365.7471 -384.1676 -375.9702
## p = 3  -337.0357 -337.0357 -340.0547 -345.3001 -367.3335 -386.7940 -378.4242
## p = 4  -322.4616 -347.1152 -347.1152 -346.5539 -367.4546 -386.1632 -377.7601
## p = 5  -319.8865 -344.6161 -342.6206 -342.6206 -365.7307 -384.5430 -376.0838
## p = 6  -341.4224 -364.7537 -363.6355 -368.6268 -368.6268 -382.5433 -374.0839
## p = 7  -340.8087 -378.3434 -377.5859 -375.6027 -376.7696 -376.7696 -375.0578
## p = 8  -354.5791 -366.2709 -364.3006 -362.7393 -375.7278 -375.7105 -375.7105
## p = 9  -350.4292 -363.9184 -362.3962 -364.7952 -372.1322 -371.7847 -378.4142
## p = 10 -368.8111 -373.3366 -372.4250 -376.3493 -377.4959 -375.8399 -375.5465
## p = 11 -344.8765 -356.5942 -354.8268 -354.7728 -365.0087 -363.9262 -365.3053
## p = 12 -374.0916 -377.3696 -379.7657 -384.3032 -383.1335 -381.9391 -381.5033
## p = 13 -370.0140 -370.4846 -373.4246 -379.0872 -379.0090 -379.7064 -383.9345
## p = 14 -390.6358 -399.7993 -404.4356 -412.6724 -410.7726 -408.9217 -407.4878
## p = 15 -399.5293 -399.7197 -409.3180 -413.4357 -413.2173 -412.9722 -413.1887
##            q = 8     q = 9    q = 10    q = 11    q = 12    q = 13    q = 14
## p = 1  -384.1193 -378.2772 -375.7594 -373.9045 -385.5903 -387.3381 -386.0167
## p = 2  -382.8545 -376.5091 -373.9810 -372.3933 -386.4067 -391.6719 -393.2430
## p = 3  -384.4024 -378.3598 -374.2608 -371.9108 -384.8451 -390.8902 -393.3453
## p = 4  -385.5532 -378.6670 -375.5295 -372.5432 -388.0965 -393.6976 -396.4510
## p = 5  -385.9886 -379.7186 -375.7508 -372.1253 -386.3250 -392.7606 -397.2054
## p = 6  -385.7800 -380.4529 -375.2765 -372.1626 -385.7889 -392.3472 -397.2233
## p = 7  -384.4661 -379.0075 -373.7403 -370.4953 -384.0889 -390.7900 -396.5567
## p = 8  -382.4798 -377.0412 -371.8235 -369.5193 -382.1253 -388.9047 -394.6737
## p = 9  -378.4142 -377.4935 -373.0680 -368.9819 -380.9882 -387.5795 -392.6889
## p = 10 -376.4238 -376.4238 -374.7640 -372.0928 -379.5501 -388.0001 -393.1039
## p = 11 -369.8459 -370.8232 -370.8232 -377.4238 -381.6288 -387.3806 -396.6590
## p = 12 -382.0895 -380.2965 -387.4949 -387.4949 -386.6923 -391.0751 -402.1894
## p = 13 -382.8989 -380.9152 -382.8796 -382.8836 -382.8836 -394.3992 -405.4274
## p = 14 -407.1036 -407.2484 -406.7846 -411.3588 -409.4708 -409.4708 -415.8596
## p = 15 -411.1888 -413.6247 -414.2691 -418.0044 -418.2690 -420.8662 -420.8662
##           q = 15
## p = 1  -406.2971
## p = 2  -407.7076
## p = 3  -405.7860
## p = 4  -404.5495
## p = 5  -402.9845
## p = 6  -403.4958
## p = 7  -401.8059
## p = 8  -400.3986
## p = 9  -398.4607
## p = 10 -398.2944
## p = 11 -402.1676
## p = 12 -405.5036
## p = 13 -406.7592
## p = 14 -422.5152
## p = 15 -421.6551
## 
## $min.Stat
## [1] -422.5152

ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, cases_growth_daily), cases_growth_daily~leaving_home_dif)

## $p
##   leaving_home_dif
## 1               14
## 
## $q
## [1] 15
## 
## $Stat.table
##           q = 1    q = 2    q = 3    q = 4    q = 5    q = 6    q = 7    q = 8
## p = 1  455.1482 433.3015 410.9211 391.1225 383.7757 325.2643 317.0773 309.5257
## p = 2  438.6280 434.9176 407.7428 391.0769 382.9822 326.1946 318.4596 310.7772
## p = 3  429.0497 429.0497 409.2341 392.5283 384.5784 327.7080 319.1513 312.0781
## p = 4  392.5450 394.4315 394.4315 391.1500 385.6264 328.5240 320.6296 312.7734
## p = 5  405.7746 405.3126 387.3259 387.3259 384.5380 330.4206 322.5490 314.4807
## p = 6  391.2806 392.5195 370.2129 369.0508 369.0508 332.4052 324.2641 315.0921
## p = 7  345.3740 342.9651 323.9536 325.5670 326.5344 326.5344 325.7480 316.7505
## p = 8  337.8587 338.6842 331.4975 333.4971 335.2338 316.6973 316.6973 318.6885
## p = 9  340.4189 341.4923 326.0076 328.0047 329.7013 314.0485 315.4048 315.4048
## p = 10 327.3648 327.7977 313.0216 312.8989 307.4921 292.7150 294.7117 296.1535
## p = 11 308.1517 308.9899 306.0758 305.7795 307.0029 299.7709 301.6611 301.4981
## p = 12 257.2831 257.5281 258.3869 259.2928 254.8132 256.3251 255.6834 256.7494
## p = 13 245.4109 242.7546 243.4022 243.8525 236.1568 237.2273 233.9605 226.8964
## p = 14 232.6548 232.9970 223.7710 220.3923 211.3452 212.4749 214.4128 214.5589
## p = 15 229.3912 225.4392 225.1676 216.3217 209.0668 208.7627 209.3096 209.1429
##           q = 9   q = 10   q = 11   q = 12   q = 13   q = 14   q = 15
## p = 1  304.6220 295.0117 289.0769 247.0309 230.1179 224.5419 205.4068
## p = 2  306.1870 296.2389 290.1364 249.0302 231.6995 225.5908 205.8936
## p = 3  306.6488 297.4570 292.0421 250.3630 233.6176 227.5879 207.7236
## p = 4  307.5799 298.2384 292.8268 251.0191 235.3207 228.6824 209.3544
## p = 5  307.5253 298.7816 294.1911 252.9023 237.3195 230.6396 211.3489
## p = 6  307.4943 297.8448 293.6651 250.4613 229.9375 226.1907 207.8137
## p = 7  309.4431 299.4357 295.6461 252.3139 230.5858 225.9718 208.4550
## p = 8  311.2740 300.7861 297.0959 253.0456 230.0256 225.7525 209.9828
## p = 9  306.5612 296.9519 294.7386 252.4602 226.8062 220.0311 209.0678
## p = 10 296.1535 297.2940 293.6219 254.4226 228.6868 221.8292 211.0547
## p = 11 296.5342 296.5342 294.1716 254.9051 230.6857 223.8239 211.0466
## p = 12 255.7004 257.0406 257.0406 254.2740 227.2089 223.0372 208.4518
## p = 13 224.5380 226.5380 225.4304 225.4304 225.9653 218.9245 206.1911
## p = 14 206.7271 207.9034 207.5801 201.3779 201.3779 203.2447 195.2065
## p = 15 210.6464 210.8963 207.1847 200.7441 195.8866 195.8866 196.5779
## 
## $min.Stat
## [1] 195.2065

testing_ardl4_growth = ardlDlm(x = scc_test_big$leaving_home_dif, y = scc_test_big$cases_growth_daily, p = 5, q = 5, remove = list(p = c(0,1,2,3)))
summary(testing_ardl4_growth)

## 
## Time series regression with "ts" data:
## Start = 6, End = 83
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.8805  -0.8568   0.0257   0.6134  12.4073 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.62782    3.16213   0.831   0.4088    
## X.4          0.25736    0.13171   1.954   0.0547 .  
## X.5         -0.16885    0.13355  -1.264   0.2103    
## Y.1          0.25162    0.10391   2.421   0.0181 *  
## Y.2          0.16684    0.10285   1.622   0.1093    
## Y.3         -0.07867    0.10026  -0.785   0.4353    
## Y.4          0.38987    0.08465   4.606 1.79e-05 ***
## Y.5          0.06010    0.08504   0.707   0.4820    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.344 on 70 degrees of freedom
## Multiple R-squared:  0.8532, Adjusted R-squared:  0.8385 
## F-statistic: 58.12 on 7 and 70 DF,  p-value: < 2.2e-16

Scatter plot to look at thresholds

scc_nolag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving_mov) %>%
  cbind(`Lag` = "0")

scc_3lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving3_mov) %>%
  cbind(`Lag` = "3")
colnames(scc_3lag)[ncol(scc_3lag)-1] = "leaving_mov"

scc_4lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving4_mov) %>%
  cbind(`Lag` = "4")
colnames(scc_4lag)[ncol(scc_4lag)-1] = "leaving_mov"

scc_5lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving5_mov) %>%
  cbind(`Lag` = "5")
colnames(scc_5lag)[ncol(scc_5lag)-1] = "leaving_mov"

scc_6lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving6_mov) %>%
  cbind(`Lag` = "6")
colnames(scc_6lag)[ncol(scc_6lag)-1] = "leaving_mov"

scc_7lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving7_mov) %>%
  cbind(`Lag` = "7")
colnames(scc_7lag)[ncol(scc_7lag)-1] = "leaving_mov"

scc_8lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving8_mov) %>%
  cbind(`Lag` = "8")
colnames(scc_8lag)[ncol(scc_8lag)-1] = "leaving_mov"

scc_9lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving9_mov) %>%
  cbind(`Lag` = "9")
colnames(scc_9lag)[ncol(scc_9lag)-1] = "leaving_mov"

scc_10lag <- scc_test_big %>%
  dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving10_mov) %>%
  cbind(`Lag` = "10")
colnames(scc_10lag)[ncol(scc_10lag)-1] = "leaving_mov"

scc_scatter <- rbind(scc_nolag, scc_3lag, scc_4lag, scc_5lag, scc_6lag, scc_7lag, scc_8lag, scc_9lag, scc_10lag)

fig_cases_growth <- 
  plot_ly (scc_scatter) %>%
    add_trace(
      x = ~leaving_mov, 
      y = ~cases_growth_daily_mov, 
      frame = ~`Lag`, 
      type = 'scatter', 
      mode = 'markers', 
      showlegend = F
    ) %>%
  animation_button(visible = F) %>%
  layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily case growth rate'), title = list(title = "SCC"))

fig_cases_growth

fig_cases_abs <- 
  plot_ly (scc_scatter) %>%
    add_trace(
      x = ~leaving_mov, 
      y = ~dcases,
      frame = ~`Lag`, 
      type = 'scatter', 
      mode = 'markers', 
      showlegend = F
    ) %>%
  animation_button(visible = F) %>%
  layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily new cases'), title = list(title = "SCC"))

fig_cases_abs

Plot the case growth rate and social distancing data

4 and 6 day lags seem the best, depending on which Granger test metric was used.

scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4_mov, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4_mov, color="Leaving home")) +
  geom_line(aes(y = dcases_mov-40, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Daily new cases, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4_mov, color="Leaving home")) +
  geom_line(aes(y = d2cases_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Change in change in cases, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

# no moving average
scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving4, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 4 days ago relative to before", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 4 Day Lag")

# 6 day
scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 6 Day Lag")

# no moving average
scc_test_big %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "Santa Clara County Case Growth and Social Distancing, 6 Day Lag")

Try to animate.

second_axis <- list(overlaying = "y",
  side = "right",
  title = "Daily case growth rate (%), 7 day moving average")

fig_animate_lines_scc <- 
  plot_ly (scc_scatter %>% mutate(date_numeric = as.numeric(date))) %>%
  add_trace(
    x = ~date_numeric,
    y = ~leaving_mov,
    frame = ~`Lag`,
    type = 'scatter',
    mode = 'marker'
    #name = "% leaving"
  ) %>%
  animation_button(visible = F) %>%
  add_trace(
      x = ~date_numeric, 
      y = ~cases_growth_daily_mov,
      frame = ~`Lag`,
      type = 'scatter', 
      mode = 'lines', 
      yaxis = "y2"
      #name = "cases"
    ) %>%
  layout(xaxis = list(title = 'date, numeric'), yaxis = list(title = 'Change in % of devices leaving home relative to before, 7 day moving average'), yaxis2 = second_axis, title = list(title = "SCC"))

fig_animate_lines_scc

San Mateo County

For now, just plotting the metric of cumulative percent positive

smc_cases_sd_daily <- bay_sd %>% 
  filter(origin_census_block_group %in% smc_blockgroups$GEOID) %>%
  group_by(date) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(
    percent_at_home = total_at_home*100/total_devices, 
    percent_leaving_home = (100 - percent_at_home),
  ) %>% 
  left_join(
    smc_cumulative_cases
  ) %>% 
  filter(date >= (min(smc_cumulative_cases$date)-10))

# get the baseline percent of people at home
pre_case_growth_at_home_smc <- bay_sd %>%
  filter(date < min(smc_cumulative_cases$date)) %>%
  filter(origin_census_block_group %in% smc_blockgroups$GEOID) %>%
  summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
  mutate(percent_at_home = total_at_home*100/total_devices, percent_leaving_home = (100 - percent_at_home))

# include change in percent leaving home
smc_cases_sd_daily <- smc_cases_sd_daily %>%
  mutate(leaving_home_dif = percent_leaving_home - pre_case_growth_at_home_smc$percent_leaving_home[1],
         log_cases = log(cumulative_cases))

smc_test_small <- smc_cases_sd_daily %>%
  mutate(
    leaving_mov = movavg(leaving_home_dif, 7, type = "s"),
    leaving4 = c(rep(NA,4), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-4)]),
    leaving4_mov = movavg(leaving4, 7, type = "s"),
    leaving5 = c(rep(NA,5), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-5)]),
    leaving5_mov = movavg(leaving5, 7, type = "s"),
    leaving6 = c(rep(NA,6), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-6)]),
    leaving6_mov = movavg(leaving6, 7, type = "s"),
    leaving7 = c(rep(NA,7), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-7)]),
    leaving7_mov = movavg(leaving7, 7, type = "s"),
    leaving8 = c(rep(NA,8), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-8)]),
    leaving8_mov = movavg(leaving8, 7, type = "s"),
    leaving9 = c(rep(NA,9), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-9)]),
    leaving9_mov = movavg(leaving9, 7, type = "s"),
    dcases = c(NA,diff(cumulative_cases)),
    leaving10 = c(rep(NA,10), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-10)]),
    leaving10_mov = movavg(leaving10, 7, type = "s"),
    dcases = c(NA,diff(cumulative_cases)),
    past_cases = c(NA, smc_cases_sd_daily$cumulative_cases[1:(nrow(smc_cases_sd_daily)-1)]), 
    cases_growth_daily = (dcases / past_cases)*100,
    cases_growth_daily_mov = movavg(cases_growth_daily, 7, type = "s")) %>% 
  filter(date >= "2020-03-01")

smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "San Mateo County Case Growth and Social Distancing, 6 Day Lag")

# moving average
smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Case Growth and Social Distancing, 6 Day Lag")

# try with percent positive cumulatively
smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6, color="Leaving home")) +
  geom_line(aes(y = perc_positive_cumulative-40, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Cumulative percent of tests that are positive")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")

smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = perc_positive_cumulative_mov-40, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Cumulative percent of tests that are positive, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")

smc_test_small %>% ggplot(
  aes(x = date)) +
    geom_line(aes(y = leaving6_mov, color="Leaving home")) +
  geom_line(aes(y = perc_positive_movavg*100-30, color = "Cases")) + 
  scale_y_continuous(sec.axis = sec_axis(~.*1-30/100, name = "Percent of tests that are positive, 7 day moving average")) + 
  scale_colour_manual(values = c("red", "blue")) +
  labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")

smc_nolag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving_mov) %>%
  cbind(`Lag` = "0")

smc_4lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving4_mov) %>%
  cbind(`Lag` = "4")
colnames(smc_4lag)[2] = "leaving_mov"

smc_5lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving5_mov) %>%
  cbind(`Lag` = "5")
colnames(smc_5lag)[2] = "leaving_mov"

smc_6lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving6_mov) %>%
  cbind(`Lag` = "6")
colnames(smc_6lag)[2] = "leaving_mov"

smc_7lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving7_mov) %>%
  cbind(`Lag` = "7")
colnames(smc_7lag)[2] = "leaving_mov"

smc_8lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving8_mov) %>%
  cbind(`Lag` = "8")
colnames(smc_8lag)[2] = "leaving_mov"

smc_9lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving9_mov) %>%
  cbind(`Lag` = "9")
colnames(smc_9lag)[2] = "leaving_mov"

smc_10lag <- smc_test_small %>%
  dplyr::select(cases_growth_daily_mov, leaving10_mov) %>%
  cbind(`Lag` = "10")
colnames(smc_10lag)[2] = "leaving_mov"

smc_scatter <- rbind(smc_nolag, smc_4lag, smc_5lag, smc_6lag, smc_7lag, smc_8lag, smc_9lag, smc_10lag)

fig_cases_growth_smc <- 
  plot_ly (smc_scatter) %>%
    add_trace(
      x = ~leaving_mov, 
      y = ~cases_growth_daily_mov, 
      frame = ~`Lag`, 
      type = 'scatter', 
      mode = 'markers', 
      showlegend = F
    ) %>%
  animation_button(visible = F) %>%
  layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily case growth rate'), title = list(title = "SMC"))

fig_cases_growth_smc